Technical Report/Research Paper No. 54 



UNITED STATES 

NAVAL POSTGRADUATE SCHOOL 




LAMINAR CONVECTIVE HEAT TRANSFER 
IN THE ENTRANCE REGION 
BETWEEN PARALLEL FLAT PLATES 



***** 

D. D. Lundberg 
and 

Dr.Jj. A. Miller 






JUNE 1965 

TECHNICAL REPORT No. 54 



literary 

u. S. Neva! Postaraduate 
Monterey, Calif oniia 



Schoo? 



LAMINAR CONVECTIVE HEAT TRANSFER 
IN THE ENTRANCE REGION BETWEEN PARALLEL FLAT PLATES 



by 



D. D. Lundberg 
and 

Dr. J. A. Miller 



U. S. Naval Postgraduate School 
Monterey, California 
19 6 5 



TM. rio, S<f 



ABSTRACT 



Heat transfer rates for laminar, convective heat transfer in the 
entrance region between parallel plates were investigated. The hydro- 
dynamic solution due to Bodia T2l was used in the solution of the energy 
equation in finite difference form on a digital computer. The thermal 
boundary conditions include: constant heat input* constant wall tempera- 
ture, one wall constant temperature and one wall insulated, and constant 
but different wall temperatures on the upper and lower walls. 

The approximate integral methods of Siegel and Sparrow T 51 , T 7 1 are 
shown to produce results that are in close agreement with the solutions 
in the present analysis for the constant heat input and constant wall 
temperature cases . 

The scope of the finite difference solution is limited to a narrow 
range of Prandtl numbers near unity, due to the small grid sizes required 
for convergence. 
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I. INTRODUCTION 



Efficient design of compact heat exchangers such as are employed for 
gas turbine regenerators, among other applications, requires maximum utili- 
zation of the high transfer rates available in the entrance region. In 
this region the velocity and temperature profiles of the entering fluid are 
undergoing rapid transition from their uniform distribution at entrance to 
the fully established distributions encountered farther downstream. The 
well known Graetz solution, in which only the temperature development is 
considered, does not yield adequate results when the velocity and tempera- 
ture distributions are simultaneously developing. This is particularly 
true if the entrance region represents a considerable percentage of over- 
all length and the two profiles develop at approximately the same rates™ 

This situation is frequently encountered with gases where the Prandtl 
number is approximately one. 

The purpose of the present analysis is to present a finite difference 
solution for the heat transfer rates between parallel flat plates for 
several different boundary conditions corresponding to various constant 
wall temperatures and constant heat input configurations. The approach 
to the problem was suggested by Miller Til, which is adapted from the num- 
erical procedures applied to the entrance region of the circular tubes by 
Kays f 2 “| * 

While an analysis such as the present one can hardly compete for 
simplicity with approximate treatments such as the Graetz analysis, or 
the several energy integral analyses commonly employed in engineering design, 
it does provide an exact solution for a limited number of cases which serves 
as a standard with which the accuracy of the various approximate methods 
may be assessed. 
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II. ANALYSIS 



A. Governing Equations 

The governing equations are the energy, momentum (Navier-Stokes) and 
continuity equations. In order to reduce the complexity and coupling of 
the equations, the following assumptions are made. 

The flow is assumed to be: 

1. steady s two-dimensional 

2. laminar 

3. incompressible 

In addition it is assumed that: 

1. thermal diffusivity (cc) is constant 

2. convective heat transfer is large compared to radiation, axial 
conduction, and viscous dissipation. 

The governing equations in reduced form may then be expressed as: 

Energy: 



u 






(i) 



Momentum: 



du , du 
u + V 

Sx ay 



- S£ + u i£u 
P 3 X 



Continuity: 



&- + 2L = o 

ax ay 



( 2 ) 



(3) 



Due to the assumption of incompressibility the above equations are 
no longer coupled. The temperature development remains dependent upon 
the velocity profile, although dependence of velocity on temperature has 
been removed. 
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The boundary conditions associated with equations (1), (2), and (3) 
for flow between parallel flat plates with uniform velocity and temperature 



at entrance are: 



1. At the entrance (x = 0) 




u = u 0 = c 2 

v = v 0 = 0 

2. At the walls (y = + d) 

v = v — 0 
w 

u = u =0 



3. On centerline 
v = v c = 0 

4. Applied thermal boundary conditions 

(1) symmetric cases (with respect to duct centerline) 
a. constant wall temperature 




b. constant heat input 



(2) asymmetric cases 

a. one wall constant temperature, one wall insulated 



w 



t » t 



w 






b. Constant wall temperature 



^ " C 6 



3 



where: f C -j 

Solution of the energy equation requires that a hydrodynamic solution be 
first obtained. 

B. Hydrodynamic Solution 

Several hydrodynamic solutions for the entrance region are available. 
Solutions have been presented by Bodoia f 3l 5 Schlichting T4]<> Schiller ^ 5l s 
and Han T61. The approximate solutions of Schiller and Han have the advantage 
of being ana lytic s and therefore, less cumbersome to use. They are, however, 
inherently less accurate than the numerical solutions presented by Bodoia and 
Schlichting. 

Bodia’s solution was obtained by evaluation of the Prandtl boundary layer 
equations in finite°dif ference form. Velocity profiles were obtained at given 
streamwise locations by utilizing matrix methods to simultaneously solve a 
column array of momentum expressions normal to the flow. 

A comparison of Schlichting 1 s series solution and the Bodoia solution 
indicates that ? in addition to the velocity gradient discontinuity in the 
Schlichting solution where the upstream and downstream solutions are joined, 
Schlichting 1 s representation of the pressure gradient, based only on the 
centerline velocity,, results in a core velocity gradient which is too large 
and a wall velocity gradient which is too small T 3 1 . 

Comparisons of the dimensionless pressure decrement versus length due 
to Bodoia^ Schlichting and Han are presented in Fig. 2. Velocity distribution 
comparisons between Bodoia and Schlichting > and Schlichting and Schiller, are 
shown in Figs. 3 and 4. 

The Bodoia solution was adopted for the present analysis. 
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C. Energy Solution 



Equations (1) and (3) may be written in dimensionless form by introduc- 
ing the following dimensionless variables^ 

T = t/t X = -\ X - = x/d U = u/u Q 

° dV R 

o e 

Y “ y/d V = v/u Q 



In dimensionless form, the equations are then: 
Energy: 



U 



+ ReviSI 

SX dY 



_ 1- 3 2 T 

PR BY^ 



Continuity: 



dx 



+ 




0 



( 4 ) 

(5) 



The following finite difference approximations are then introduced: 

ST _ AT Tx+l s y - Tx,y 
SX ~ AX 2Tx 



ST _ At _ Tx,y+1 = Tx,y=l 
a Y Ay " o-A^ 



a^T = Tx,y+1 ~ 2Tx,y + Tx,y=l 

ax a 



SU _ Au _ Ux+1 ,y-Ux-l ,y 
SX Ax ~ 2 A X 



SV _ Av _ Vx,y=Vx,y-l 

sy ay — 2Ty (10) 



Note that (6) is a forward difference and (10) is a backward difference, 
whereas the remaining equations are of the central defference type. These 
have been expressed in this manner for simplicity and convenience. The larger 
the truncation error normally associated with the forward or backward difference 
is not significant here since practical grid dimensions require that: 
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Ay » Ax 



and, due to the nature of the flow; 

Au Ay 
Ax >:> AY 

Introducing (6) through (10) into (4) and (5) and solving for the ap~ 
propriate term yields: 



Tx+1 ,y = 



L uprSy 2 + 

U2AY J 



2AX 

L 1 upray 2 J Tx >y +1 + 



I 

L upray^ 



AXReV 

U2AY 



j Tx,y+1 



(ID 



v x,y - Vx,y-1 + Re2 Ax L Ux_1>y ” u x+l,y j (12) 

which contribute the computing equations. 

Note that the Reynolds number (R e ) s although appearing in the above 
equations, remains only for convenience and could be eliminated. It is no 
longer a parameter in the set of equations haveing been included in the 
dimensionless length variable (X) . 

A mixed mean temperature (t m ) may be expressed in the following way: 



m 



j CpputdA 
j^CpPudA 



(13) 



Removing the constants from within the integral s and simplifying, a 
dimensionless mixed mean temperature may be written in finite difference 
form: 
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Ux,y Tx,y 



y=l 

l 



m y=o 



y-i 

I Ux,y 



( 14 ) 



y» 0 

Due to the no°slip condition imposed upon viscous flows, heat trans 
fer to or from the fluid at the wall may be expressed in terms of con- 
duction through a thin film of stagnant fluid. Thus: 



q" * q/A - « (f^ 



(15) 



which must equal the heat flux due to convection, 

q" = q/ A - h < t w - tin) 

Defining the local Nusselt: number as: 



Combining (14) and (15), 



N = M 
u x k 



N n „ = Qy)y* 



(16) 



(17) 



Ux 



(18) 



fc w " fc m 



A Nusselt number based on the hydraulic diameter (d) instead of half 
the gap width (d) in finite difference form becomes: 



N, 



u x 



" T w-Ay^ 

AY 

T - T 
w m 



(19) 



The hydraulic diameter was introduced in (19) to enable comparisons 
with previous solutions. The conversion is: 

D = 4d 

In addition, for the asymmetric constant wall temperature case, it 
is convenient to use a dimensionless parameter other than the Nusselt 



7 



number. A dimensionless heat flux parameter may be defined by using (15) 
above: 




D. Thermal Boundary Conditions in Finite Difference Form Constant Wall 
Temperature (Symmetric Case ). Since the temperature distribution 

is symmetrical with respect to the centerline, computations were performed 
in only the upper half of the duct. A coordinate system was adapted with 
the X axis on the centerline of the channel. The thermal boundary condi- 
tion at the wall was applied as: 



t /t 
w o 



T 

w 



* N 



where N is a positive integer. 

Constant Heat Input (Symmetric Case) 0 With the axis system as above, 
the slope of the thermal profile was maintained constant at the wall by 
adding a constant to the calculated temperature one step from the wall, 
thus : 



where : 



T = T ~ A 
w w Ay 



+ M 



M = (AT/AY) w AY 

Constant Wall Temperature with One Wall Insulated . For this case the 
coordinate system was redefined so the X axis coincided with the lower wall. 
The lower wall boundary condition was taken to be: 



t / 1 
W]^ o 



T = N 
W 1 



while zero slope at the upper wall was obtained by equating the wall 
temperature to the calculated temperature one step from the wall: 
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tw 2 y/ft o ~ T w2 T w2 



Ay 



Constant Wall Temperature (Asymmetric Case) . With the axis system 
the same as for the insulated wall case, thermal boundary conditions were 
established as: 



^ 1^0 ~ T «1 ~ N 1 



*^w2 *** ^2 



where N-^ and Ng are positive integers and N-^ N 2 « 
Noting that: 







V^ 




1 



1 



Values of - I/N 2 - 1 of 1/2, 1/3, 1/4 were investigated. 

E . Computational Procedure 

The finite difference solution was obtained by using a CDC 1604 
digital computer to solve the equations at intervals of AX - 0.001 and 
AY * 0.1. Velocity distributions were obtained by interpolating plot- 
ted curves of Bodoia f s tabulated results. 

Since the entire solution is dependent upon the first calculations 
at the entrance, starting values were obtained by reducing the grid size 
to AX = 0.0001 and Ay * 0.05 for the first 20 streamwise stations from 
X a 0 to 0.0020. 

The continuity equation (2) was evaluated and the cross-wise velocity 
component (V) was introduced into the energy equation (11) at each grid 
point. The energy solution was then obtained at intervals of Ay between 
centerline and the wall at each streamwise station. 
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Solutions for each thermal boundary condition were obtained for 
Prandtl numbers of: 0,5, 0.7, 1.0, lo6, 3.2 and 10. 

Solutions for Prandtl numbers less than 0.5 are excluded by the grid 
size. This may be seen from inspection of the coefficient of the second 
temperature term in (11), which indicates the interdependence of PR, Ax, 
and AY. In order to obtain convergence, PR must be large enough to pre- 
vent this term from becoming negative. Thus low Prandtl numbers, ap- 
proaching those of liquid metals, cannot be effectively handled by finite- 
difference methods owing to the large number of calculations required by 
the small grid. 

The computer program, written in Fortran 60, is included in the 
Appendix. 
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Ill RESULTS 



Tables I through VI contain computed and extrapolated values of the 
heat transfer rates as a function of distance from the entrance for each 
of the thermal boundary conditions for a Prandtl number of 0.7. Figures 
5 through 29 contain curves of the transfer rate versus downstream posi- 
tion and temperature and velocity profiles for Prandtl numbers of 0.5, 

0.7, 1.0, 1.6, 3.2 and 10. 

With the exception of the Asymmetric Constant Wall Temperature Case, 
which is expressed in terms of the dimensionless heat flux parameter, q*; 
the heat transfer rates are presented as local Nusselt number variation 
with the length parameter, L. Whereas, the computations were carried out 
using the dimensionless length parameter X, which is normally employed 
for hydrodynamic developments; L, based on the hydraulic diameter is 
customarily employed in heat transfer. The conversion is: 

X = 16L 

The use of all finite starting length in the finite difference equa- 
tions results in a finite initial value for the Nusselt number at X - 0 
rather than the theoretically infinite value at the entrance. This causes 
a nearly flat slope near the entrance and an artificial inflection point 
in the Nusselt number versus length curve. The inflection point occurs 
progressively further downstream with the increasing Prandtl number, 
thereby limiting the useful span of the curves. A reduction of grid size 
from = .001, Ay = .1 to AX = .0001, AY = .05 was observed to double 
the initial Nusselt number from 40 to 80, and move the inflection point 
upstream. This increased the useful span. However, the approach obviously 
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reaches a practical limit very rapidly due to the number of calculations 
and the accuracy required of the hydrodynamic solution very close to the 
entrance. 

In general , for the three thermal boundary conditions in which the 
Nusselt number was used as a heat transfer parameter, the curves for 

Prandtl numbers of 0.5, 0.7, 1.0 and 1.6 behave as expected for values of 

-A — 3 

L >10” . For Prandtl numbers of 3.2 and 10, L * 10 is the approximate 

lower limit of validity. 

Constant Heat Input . A comparison of local Nusselt numbers obtained 
in the present analysis with those reported by Siegel and Sparrow T7l, and 
Han T61, contained in Fig. 5, indicate close agreement with the results of 
Siegel and Sparrow. Their solution employed a simplified energy equation 
in which the temperature distribution in the boundary layer was expressed 
as a series of polynomials in the transverse coordinate using the down- 
stream station as a parameter. The hydrodynamic solution employed 
Schilleir 1 s approximation. 

The comparison with Han's integro-numerical solution is not as close 
as might appear in Fig. 5 since this curve is for a Prandtl number of 0.8 
and should lie above the other two curves which are for a Prandtl number 
of 0.7. 

The computed curves of local Nusselt number versus length contained 

in Fig. 6, for Prandtl numbers of 0.5, 0.7, 1.0, 1.6, 3.2 and 10, may 

-3 

be approximated within 10 percent, for 0.5 ^PR^1.6, and L >10 by the 
following: 



N u = 8.24 + 




x 



1 + .0178 (PR/L)-6 
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Velocity and temperature profiles for this boundary condition are con- 
tained in Figs. 13 through 16. 

One Wall faonstant Temperature - One Wall Insulated . The results of this 

computation are presented in Fig. 17. The local Nusselt number may be approx- 

-3 

imated within 10 percent for 0.5sPR^1.6, and L >10 by: 



N u = 4.84 + 
u x 



1 + 



PR 

. 0155( /L) 

. 012 ( PR /L) 



The velocity and temperature distributions are illustrated in Figs. 

18 through 21. 

Constant Wall Temperature (Asymmetric Case) . Dimensionless heat flux 
values for this boundary condition are presented in Figs. 22, 23 and 24 
for ^1 ^ ratios of 1/2, 1/3, and 1/4. 



Note that the upper (hotter) wall curve is smoothly asymptotic to the 
value predicted for the steady state condition. The lower curve exhibits a 
mush shallower initial slope, passes through an inflection point , and then s 
approaches a negative asymptotic value of the same magnitude and at the 
same rate as the upper curve. 

Velocity and temperature profiles for N n ^l are presented in 

~~r = 1/4 

Figs. 25 through 28. 2 

The behavior of the local Nusselt number for this asymmetrical case 
is depicted in Fig. 29. The lower curve, representing the local Nusselt 
number variation for the lower (cooler) wall, does not depict the actual 
physical situation. 
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Both the denominator and the numerator of; 



4(T 



w 






N 



ay 



T - T 
w m 



( 18 ) 



approach zero at different rates. The curve tends to positive infinity 
as the denominator approaches zero, returns from negative infinity, and 
becomes zero as the numerator goes to zero; and finally approaches the 
actual steady asymoptotic value. 
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IV. CONCLUSIONS 



From the results of this study it may be concluded that: 

1. The finite difference method, as employed in this analysis, yields re- 
sults which are consistent with previous solutions for other geometries T2l 
and which are asymptotic to the known fully developed values* The 
method is, however, limited to a small range of Prandtl numbers near 
unity by the small grid sizes necessary outside this range* For small 
Prandtl numbers the grid must be very small to obtain convergence, 

while for large Prandtl numbers practical grid dimensions result in 
significant inaccuracies near the entrance due to the "finite starting 
length 11 • 

2. The results of the approximate integral methods suggested by Siegel and 
Sparrow T7] for constant heat input, and by Sparrow T 5 1 for constant 
wall temperature, compare very favorably with the finite difference 
solutions reported here. This additional substantiation of the methods 
suggested by these authors is significant because of the practical 
utility of their approach, which requires relatively few calculations, 
is capable of handling a wide range of Prandtl numbers, and can pre- 
dict heat transfer rates very close to the entrance* 
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TABLE I 



CONSTANT HEAT INPUT 
PR = 0.7 



X 


L 


V 


T 

w 


.0016 


.0001 


46.2 


1.09 


.0032 


.0002 


30.0 


1.14 


.0064 


.0004 


23.7 


1.18 


.0128 


.0008 


18.0 


1.24 


.0160 


.0010 


16.6 


1.26 


.0320 


.0020 


12.9 


1.35 


.0640 


.0040 


10.3 


1.47 


.1280 


.0080 


8.87 


1.62 


.1600 


.0100 


8.62 


1.67 






Asymptotic to 8.24 
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TABLE II 





CONSTANT WALL TEMPERATURE 






(SYMMETRIC CASE) 






PR = 0.7 




X 


L 


% 


.0016 


.0001 


34. 


.0032 


.0002 


26. 


.0064 


.0004 


19. 


.0128 


.0008 


15. 


.0160 


.0010 


13. 


.0320 


.0020 


10. 


.0640 


.0040 


8. 


.1280 


.0080 


7. 


.1600 


.0010 


7. 



Asymptotic to 7. 



0 

0 

9 

3 

7 

5 

58 

78 

63 

60 



18 



TABLE III 





ONE WALL CONSTANT TEMPERATURE - ONE 


WALL INSULATED 






PR = 0.7 






X 


L 




T 

w 


.0016 


.0001 


33.0 


1.000 


.0032 


.0002 


24.0 


1.000 


.0064 


.0004 


19.2 


1.000 


.0128 


.0008 


13.0 


1.000 


.0160 


.0010 


11.6 


1.000 


.0320 


.0020 


9.05 


1.000 


.0640 


.0040 


7.20 


1.000 


.1280 


.0080 


5.93 


1.001 


.1600 


.0100 


5.55 


1.021 


.4000 


.0250 


4.96 


1.139 




Asymptotic to 


4.84 


2.000 
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TABLE IV 



CONSTANT WALL TEMPERATURE 
(ASYMMETRIC CASE) 



N x - 1 - 1/2 
N 2 - 1 



X L 



0032 


.0002 


0064 


.0004 


0128 


.0008 


0160 


.0010 


0320 


.0020 


0460 


.0040 


1280 


.0080 


1600 


.0100 


3200 


.0200 


6400 


.0400 


9600 


.0600 



Asymp 



PR = 0.7 



* * 



91 


q 2 


26.2 


57.5 


17.3 


41.4 


11.9 


28.7 


10.7 


25.9 


8.08 


18.7 


6.18 


13.7 


4.73 


10.1 


4.24 


9.15 


2.24 


6.35 


-.19 


3.81 


-1.23 


2.76 


-2.00 


+2.00 
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TABLE V 



X 

.0032 

.0064 

.0128 

.0160 

.0320 

.0640 

.1280 

.1600 

.3200 

.6400 

.9600 



CONSTANT WALL TEMPERATURE 
(ASYMMETRIC CASE) 



N L - 1 = 1/3 PR = 0. 




L 

.0002 

.0004 

.0008 

.0010 

.0020 

.0040 

.0080 

.0100 

.0200 

.0400 

.0600 

Asymptotic to 



* * 





q 2 


28.0 


86.0 


17.3 


62.0 


11.8 


43.0 


10.8 


38.9 


8.09 


28.1 


6.18 


20.6 


4.69 


15.2 


4.15 


13.8 


1.65 


9.84 


- 1.58 


6.41 


- 2.97 


5.02 


- 4.00 


+ 4.00 
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TABLE VI 





CONSTANT WALL TEMPERATURE 








(ASYMMETRIC CASE) 








N x - 1 = 1/4 PR = 0.7 








i— 1 
1 

<N 

S3 


* 


* 


X 


L 


91 


*2 


.0032 


.0002 


28.2 


116. 


.0064 


.0004 


16.8 


82. 


.0128 


.0008 


12.2 


57. 


.0160 


.0010 


10.5 


50. 


.0320 


.0020 


8.09 


37. 


.0640 


.0040 


6.17 


27. 


.1280 


.0080 


4.66 


20. 


.1600 


.0100 


4.06 


18. 


.3200 


.0200 


1.06 


13. 


.6400 


.0400 


-2.96 


9. 


.9600 


.0600 


-4.71 


7. 




Asymptotic to 


-6.00 


+6. 



0 

9 

7 

4 

5 

5 

3 

4 

3 

03 

27 

00 
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FIGURE 6 



TEMPERATURE AND VELOCITY PROFILES 
X= 0.Q50 i L = 0.003125 
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FIGURE 9 
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FIGURE 13 
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FIGURE 14 
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FIGURE 15 
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TEMPERATURE AND VELOCITY PROFILES 
ONE WALL CONSTANT TEMPERATURE 
ONE WALL INSULATED 
X= 0.005 L = 0.0003125 
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FIGURE 19 
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, ONE WALL INSULATED 
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FIGURE 20 
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FIGURE 21 
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FIGURE 26 

TEMPERATURE AND VELOCITY PROFILES 
CONSTANT WALL TEMPERATURE 
(ASYMMETRIC CASE) 

X= 0.050 L =0.003125 
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FIGURE 28 

TEMPERATURE AND VELOCITY PROFILES 
CONSTANT WALL TEMPERATURE 
(ASYMMETRIC CASE) 

X* 1.00' L= 0.0625 
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FIGURE 29 



* LOCAL NUSSELT NUMBER 
CONSTANT WALL TEMPERATURE 
(ASYMMETRIC CASE) 

N, -I _ I 
N 2 -I ” 4 

PR = 0.7 
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hnnnnnnrinnnnn 



. . JOB0191F.LUNDBERG 

PROGRAM HEATX23. 
i CONSTANT HEAT INPUT 
BODOIA VELOCITY DATA 

THIS PROGRAM USES A REDUCED GRID OF. X = .000 1 , Y = . 05 
FOR THE FIRST 20 STREAMWISE CALCULATIONS. 

GRAPH OUTPUTS OF NUSSELT NUMBER VERSUS LENGTH PARAMETER 
FOR' PRANDTL NUMBERS OF 0 . 5 ,0 . 7 » 1 . 0 » 1 . 6 , 3 . 2 , AND 10. • 

ARE OBTAINED TO 50 PER CENT OF THE HYDRODYNAMIC , - 

DEVELOPMENT. PRINT OUTPUT INCLUDES NUSSELT NUMBER. 

SOPE OF THE TEMPERATURE PROFILE AT THE WALL , AND THE 
SLOPE OF THE TEMPERATURE PROFILE AT THE WALL, AND THE MEAN 
MIXED MEAN TEMPERATURE. IN ADDITION, VELOCITY AND 
TEMPERATURE VALUES ARE PRINTED AT EACH GRID POINT AT 
SELECTED STREAMWISE LOCATIONS. 

DIMENSION TOO, 23) ,U( 30,23 ) ,V< 30,2 3) »Y 123) ,X I 500) 
1»XN(500)»ITITLE(12) 

READ 500, ( I TITLE! I ) , 1=1.6) 

READ 500 » ( I T I,TLE( I ) » 1=7,12) 

500 FORMAT (6A8) 

DO 750 KK= 1,6 

GO T0(40, 41*42,43. 44*45), KK N " 

40 PR =0.7 
GO TO 46 

41 PR=0 • 5 
GO TO 888 

GO TO 46 • 

42 PR= 1 • 0 
GO TO 46 

43 PR = 1 .6 
GO TO 46 

44 PR = 3 .2 
GO TO 46 

45 PR= 1 0. 

46 PRINT 523 
523 FORMAT (1H1) 

PRINT 17, 'PR 

17 FORMAT (///// 10X * 17H PRANDTL NUMBER =,F5.3) * 

L= 1 ■ 

GO TO 150 

700 READ 701, (<U(I,J), J=13,22). 1=2,22) 

701 FORMAT (10F6.4) 

DO 710 1=2,22 
DO 711 J= 1 , 12 

711 U( I , J) =U( I *1-3) ' 

710 CONTINUE 
DELX=. 000 1 
DELY=0 .05 
R=1000* 

DO 110 L=2 ,21 
150 UTM=0.0 

um=o.o . 

XL = L 

X(L)=(XL-1.0)*0.0001 
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1 



DO 111 J = 3 * 22 
XJ = J 

Y(,L) = (XJ-3. 01*0.05 
Y( 231=1.0 
U( l.Jl-l. 

U( 1»23)=1.0 

V(1,J-1)=0.0 . ' 

V(1»J)=0. 

T(1,J)=1.0 ‘ 

T(1,J+1 1=1.0 
T<1»J-1)=1.0 
T ( L » 1 ) =T( L» 5 ) 

. U ( L » 1 ) =U ( L » 5 ) 

U(L, 231=0.0 
U(L,2)=U(L,4) 

V( L » 3 ) =0. 0 
V(L»2)=V(L»4) 

C THERMAL BOUNDARY CONDITIONS 

T ( L » 23 ) = T(L»2,2)+.05 
T ( L » 2 ) =T ( L » 4 ) 

IF(L-l) 109,109,108 

108 V ( L , J 1 = V<L,J-1) + DELY/(2.*DELX*R) * ( U < L-l , J 1 -U ( L + l , J ) > 

109 A= DELX/ (U(L,J)*PR*DELY**2) 

V(L, 231=0.0 

B= ( DELX*R*V ( L , J ) ) / ( 2 .*DELY*U ( L , J 1 ) 

. T(L+1,J) = (A-B)*T.<L,J + 1 ) + ( l.-2.*Al*T(L,J) + (A+B)*T(L,J-l s > 

I F ( J-3 1 113,113,112 - 

113 UTM=UTM+»5*U(L,J)*T ( L,J) 

UM = UM+, 5*U ( II » J 1 

GO TO 114 

112 UTM=UTM+U(L,J)*T(L,J) 

UM=UM+U(L,J) 

114 TM=UTM/UM 
111 CONTINUE 

Tl= ( T ( L * 23 ) -T ( L» 22 ) 1 / DELY 
T2= T ( L , 23 1 - TM 
XN(L)=4.0*T1/T2 , 

131 PRINT 118, X(L) ♦ X N ( L 1 »T2»T1,TM 

118 FORMAT (// 4H X = ♦ F6 . 4 , 5X ♦ 14H NUSSELT NO. =,F8.5»5X,8H TW-TM =, 
1F8.5,5X,8H SLOPE =,F8.5,5X,5H TM =,F8.5) 

PR I NT 1 19 ♦ ( T ( L , J 1 ♦ J = 1 , 1 2 ) 

PR I NT 1 19 ♦ ( T ( L ♦ J 1 , J=13 ,23 1 ' ' . 

119 FORMAT </9H TEMP. =,12F7.4) 

PRINT120,- ( U ( L » J ) , J=1 » 1 2 ) 

PRINT120, ( U ( L , J ) , J = 1 3 , 23 ) 

120 FORMAT (/ 9H X VEL. =,12F7.4) 

PRINT121, ( V ( L , J ) ♦ J=1 , 12 ) 

PRINT121, ( V ( L , J ) , J = 1 3 , 23 ) 

121 FORMAT (/ 9H Y VELv =,12p7.4) 

IF(L-l) 700,700,110 

110 CONTINUE 

CALL DRAW ( 2 1 ♦ X ,XN ,MOD » 0 , LAB » I T I TLE , . 04 ♦ 10 . , 0 , 0 , 0 ,0 , 6 ♦ 6*1, LAST) 
DO 720 N=2 ,11 
M=2*N~ 1 
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U( 1 »N)=U( 1 1 » M ) 

U ( 2 » N ) =U< 2 1 »M ) 
T(2*N)=T(21*M) 

. . PRINT 771, M, U(2»N) ,T(2»N) 
771 FORMAT ( I10#2F10.*3 ) 

720 CONTINUE 
85 L = 2 

DO 10 M = 2 ,250 
70 XM=M 
Y=0 • 0 
DELY=0.1 
DELX=0 • 00 1 
X ( M ) =DELX*XM 
I = L + 1 
' N3N=M 
N = L + 1 



82 I F ( M- 99) 72,73,73 

72 READ 1 (U( I ,11) ,U( I, 10 ) *U( 1,9) ,U(I ,8) »U( I ,7) ,U( I ,6) 
1 , U ( 1,5) »U( I ,4) ,U( 1,3) ,U( 1,2) ) 

1 FORMAT (10F6.4) 

73 IFIM-99 ) 50,51,52 
52 I F ( M-l 49 ) 51,53,54 
54 I F ( M-199 ) 53,55,55 



51 



UU=XM- 
U ( N , 2 ) 

U ( N , 3 ) 

U ( N , 4 ) 

U ( N , 5 ) 

U ( N , 6 ) 
U< N » 7 ) 

U ( N , 8 ) 

U ( N , 9 ) 

U ( N , 10 ) 
U ( N , 1 1 ) 



99. 

■=1.4388 + 
=1.4292+ 
=1.3993+ 
=1.3454+ 
=1.2628- 
=1.1467- 
= .9938- 
= .8022- 
= .5720- 
= .3042- 



•000740*UU 

• 000674*UU 
•000492*UU 

• 0002 38*UU 
•000034*UU 
• 000262*UU 
.000410*UU 
•000452*UU 
•000388*UU 

• 0002 30*UU 



GO TO 50 
53 UU=XM-149. 

. U ( N , 2 ) =1.4758+.000290*UU 
U ( N , 3 ) =1.462 9+ .000266 *UU 
U ( N , 4 ) = 1 ,4239+ .000194 *UU 
U ( N , 5 ) =1.3573+.000092#UU 
U ( N ♦ 6 ) =1.2611-. 000014*UU 
U ( N , 7 ) =1. 1336-.000104*UU 
U ( N , 8 ) = ,9733-.000160*UU 
U ( N , 9 ) = . 7796- .0001 76*UU 
U ( N , 10 ) = . 5526-.000150+UU 
U ( N , 1 1 ) = . 2926-.000092*UU 



GO TO 50 
55 UU=XM-199. 

U(.N,2) =1.4903+. 0000121*UU 
U ( N , 3 ) =1.4762+.0000110*UU 
U( N ,4) =1.4336+. 0000 080*UU 
U ( N * 5 ) =1.3619+. 0000026*UU 
U ( N , 6 ) =1.2604-. 0000005*UU 
U ( N , 7 ) =1. 1284-.0000042*UU 
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U ( N » 3 ) = • 9 65 3- • 0000066*UU . 

U ( N » 9 ) = .7708-.0000072*UU 
U (, N » 1 0 ) = .5451-.0000063*UU 
U(N,11)= .2880-.0000475*UU 
50 UTM=0.0 
UM=0.0 

DO 11 J = 2 » 1 1 ' . 

X J= J 

V( J)»DELY*XJ-.2 • 7 

U< L , 12 ) =0.0 
V( L » 2 ) =0 .0 
V(L,1)=V(I»3) 

Y ( 12 ) = 1 .0 

T(L»12)=T(L»11)+0.1 

T(L»1)=T(L»3) 

U(L*1)=U(L*3) 

8 V ( L » J ) = V(L.J-l) + DELY/(2.*DELX*R) * ( U( L-l » J )-U ( L + l . J ) ) 

_A = DELX/ <U( L»J)*PR*DELY**2 ) 

, V ( L * 12 ) =0 • 0 , 

8= ( DELX*R*V ( L * J ) ) / ( 2 .*DELY*U ( L » J ) ) 

TtL+1, J)= ( A-8)*T(L,J+1 ) + ( l.-2.*A)*T<L,J > + ( A+B)*T(L, J-l ) 

I F ( J-2 ) 13.13*12 

13 UTM = UTM+. 5*U ( L » J ) *T ( L » J ) 

UM=UM+.5*U( L. J) 

GO TO 14 

12. UTM=UTM+U( L»J)*T(L»J) • , > 

UM=UM+U ( L » J ) X, 

14 TM=UTM/UM 
11 CONTINUE 

Tl= ( T(L»12)-T(L»11) ) / DELY 
T2= T ( L » 12 ) -TM 
XN(M)=4.0*T1/T2 
I F ( N3N-50 ) 31.31,30 

30 IF( (N3N/25)*25-N3N) 15.31.31 

31 PRINT 1 8 ♦ X < M ) »XN(M) .T2.T1.TM 

18 FORMAT (// 4H X = » F6 . 3 »'5X , 14H NUSSELT NO. =,F8.5»5X»8H TW-TM =, 
1F8.5.5X.8H SLOPE =.F8.5,5X,5H TM =,F8.5) 

PRINT 19, (T (L.J) , J = 1 ,12) 

FORMAT ( /9H TEMP. =,12F7.4) 

PRINT 20, ( U ( L » J ) ♦ J = 1 » 1 2 ) 

FORMAT (/ 9H X VEL. =»12F7.4) 

PRINT 21, ( V ( L » J ) ♦ J = 1 ♦ 1 2 ) 

FORMAT ( / • 9H Y VEL. =,12F7.4) 

DO 60 J = 2 » 11 
T(L,J)=T(L+1»J) 

U( L-l, J)=U(L,J) 

U( L » J ) =U( L+l , J) 

CONTINUE '■ 

I F ( KK-6 ) 760,761,761 
MOD = 3 
GO TO 762 
MOD = 2 

CALL' DRAW (497.X, XN, MOD, 0, LAB, I TITLE,. 04, 10., 0,0, 0,0,6, 6,1 .LAST ). 
CONTINUE 
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